Import shp files, remove polygons less than 1ha
library(tidyverse)
library(sf)
library(future.apply)
# Import reef polygons, create union (or returns multiple intersections
geomorphic_seeding <- st_read("/Users/rof011/spatialtools/apps/remove-zones/www/geomorphic.geojson") %>%
mutate(area=as.numeric(st_area(.))) |>
st_transform(20353) |>
filter(area > 10000) |>
mutate(n=round(as.numeric(area)/1000)) |>
st_transform(4326) |>
st_crop(st_bbox(c(xmin = 145.42, xmax = 145.48, ymin = -14.72, ymax = -14.64))) |>
st_transform(20353)
## Reading layer `Lizard_Eyrie' from data source
## `/Users/rof011/spatialtools/apps/remove-zones/www/geomorphic.geojson'
## using driver `GeoJSON'
## Simple feature collection with 736 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 145.3519 ymin: -14.73534 xmax: 145.4974 ymax: -14.64425
## Geodetic CRS: WGS 84
library(tmap)
habitat_pal <- c(
"Plateau" = "cornsilk2",
"Back Reef Slope" = "darkcyan",
"Reef Slope" = "darkseagreen4",
"Sheltered Reef Slope" = "darkslategrey",
"Inner Reef Flat" = "darkgoldenrod4",
"Outer Reef Flat" = "darkgoldenrod2",
"Reef Crest" = "coral3",
"Shallow Lagoon" = "turquoise3",
"Deep Lagoon" = "turquoise4"
)
tmap_mode("plot")
#tm_basemap("Esri.WorldImagery", alpha=0.6) +
tm_shape(geomorphic_seeding) +
tm_polygons(fill="class",
fill_alpha=0.5,
fill.scale=tm_scale_categorical("habitat_pal"))
Set seed points
set.seed(101)
#remotes::install_github("r-tmap/tmap.deckgl")
# lapply sampling process
seeded_points <- st_sf(do.call(rbind, lapply(seq_len(nrow(geomorphic_seeding)), function(i) {
polygon <- geomorphic_seeding[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
}))) |> st_transform(20353)
tm_basemap("Esri.WorldImagery") +
tm_shape(seeded_points) +
tm_dots(size=0.01,
fill="white",
fill_alpha=0.5)
# Define the width and length of the rectangle
width <- 100 # Example width (in the same units as your CRS)
length <- 100 # Example length (in the same units as your CRS)
# Use lapply to create the buffered polygons
buffered_polygons <- lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_sf <- sf::st_sfc(buffered_polygons, crs = sf::st_crs(seeded_points)) |> st_as_sf()
inset <- st_bbox(c(xmin = 145.44, xmax = 145.455, ymin = -14.697, ymax = -14.692)) |>
# st_as_sfc() |>
# st_as_sf() |>
# st_set_crs(20353) |>
st_bbox()
tm_basemap("Esri.WorldImagery") +
tm_shape(buffered_polygons_sf,
bbox=inset,
crs=20353) +
tm_polygons(fill=NA,
col="white",
lwd=0.8)
calculate intersections
# Intersect the plots with the reef polygons, calculate overlap
intersections <- st_intersection(buffered_polygons_sf, st_union(geomorphic_seeding)) %>%
mutate(intersection_area = st_area(.)) %>%
mutate(percentage_overlap = (as.numeric(intersection_area)/10000) * 100)
ggplot() + theme_bw() +
geom_histogram(data=as.data.frame(intersections), aes(x=percentage_overlap), alpha=0.2, fill="turquoise", color="black", linewidth=0.2)
view in space
library(tmap)
tmap_mode("plot")
tm_basemap("Esri.WorldImagery") +
tm_shape(intersections |> filter(!percentage_overlap > 99),
bbox=inset) +
tm_polygons(fill=NA,
lwd=0.5,
col_alpha=0.5,
col="darkred") +
tm_shape(intersections |> filter(percentage_overlap > 99),
bbox=inset) +
tm_polygons(fill=NA,
lwd=1,
col="white")
time code across ~8000 points
library(tictoc)
tic()
# lapply sampling process
seeded_points <- st_sf(do.call(rbind, lapply(seq_len(nrow(geomorphic_seeding)), function(i) {
polygon <- geomorphic_seeding[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
}))) |> st_transform(20353)
# Define the width and length of the rectangle
width <- 100 # Example width (in the same units as your CRS)
length <- 100 # Example length (in the same units as your CRS)
# Use lapply to create the buffered polygons
buffered_polygons <- lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_sf <- sf::st_sfc(buffered_polygons, crs = sf::st_crs(seeded_points)) |> st_as_sf()
# Intersect the plots with the reef polygons, calculate overlap
intersections <- st_intersection(buffered_polygons_sf, st_union(geomorphic_seeding)) %>%
mutate(intersection_area = st_area(.)) %>%
mutate(percentage_overlap = (as.numeric(intersection_area)/10000) * 100)
toc()
## 12.138 sec elapsed
parallel
options(future.rng.onMisuse = "ignore")
library(future.apply)
tic()
# Set up parallel processing plan
plan(multisession, workers = availableCores())
# Parallel lapply sampling process using future_lapply
seeded_points_list <- future_lapply(future.seed=NULL, seq_len(nrow(geomorphic_seeding)), function(i) {
polygon <- geomorphic_seeding[i, ]
points <- st_sample(polygon, size = polygon$n, type = "random")
st_sf(class = polygon$class, geometry = points)
})
# Combine the list into a single sf object
seeded_points <- do.call(rbind, seeded_points_list) |> st_sf() |> st_transform(20353)
# Define the width and length of the rectangle
width <- 100 # Example width (in the same units as your CRS)
length <- 100 # Example length (in the same units as your CRS)
# Use future_lapply to create the buffered polygons
buffered_polygons_list <- future_lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
return(polygon)
})
# Combine the results into a single sf object
buffered_polygons_sf <- sf::st_sfc(buffered_polygons_list, crs = sf::st_crs(seeded_points)) |> st_as_sf()
# Intersect the plots with the reef polygons, calculate overlap
intersections <- st_intersection(buffered_polygons_sf, st_union(geomorphic_seeding)) %>%
mutate(intersection_area = st_area(.)) %>%
mutate(percentage_overlap = (as.numeric(intersection_area) / 10000) * 100)
# Reset the plan to sequential processing after completion (optional)
plan(sequential)
toc()
## 28.88 sec elapsed
now rotate polygons and loop the code
# library(sf)
# library(spatialEco)
# library(future.apply)
#
# # Set up parallel processing plan
# plan(multisession, workers = availableCores())
#
# tic()
#
# # Define rotation angles
# angles <- seq(0, 355, 5)
#
# # Parallelize the outer loop using future_lapply
# results_list <- future_lapply(angles, function(angle) {
#
# # Sequential sampling process within the loop
# seeded_points_list <- lapply(seq_len(nrow(geomorphic_seeding)), function(i) {
# polygon <- geomorphic_seeding[i, ]
# points <- st_sample(polygon, size = polygon$n, type = "random")
# st_sf(class = polygon$class, geometry = points)
# })
#
# # Combine the list into a single sf object
# seeded_points <- do.call(rbind, seeded_points_list) |> st_sf() |> st_transform(20353)
#
# # Define the width and length of the rectangle
# width <- 100 # Example width (in the same units as your CRS)
# length <- 100 # Example length (in the same units as your CRS)
#
# # Use lapply to create the buffered polygons
# buffered_polygons_list <- lapply(seq_len(nrow(seeded_points)), function(i) {
#
# # Extract the coordinates of the current point
# x <- sf::st_coordinates(seeded_points)[i, 1]
# y <- sf::st_coordinates(seeded_points)[i, 2]
#
# # Set parameters for the rectangle around the point
# x_min <- x - (width / 2)
# x_max <- x + (width / 2)
# y_min <- y - (length / 2)
# y_max <- y + (length / 2)
#
# # Create the rectangular polygon
# polygon <- sf::st_polygon(list(rbind(
# c(x_min, y_min),
# c(x_min, y_max),
# c(x_max, y_max),
# c(x_max, y_min),
# c(x_min, y_min)
# )))
#
# # Convert to sf object for rotation
# polygon_sf <- st_sf(geometry = st_sfc(polygon, crs = st_crs(seeded_points)))
#
# # Rotate the polygon by the specified angle using spatialEco::rotate.polygon
# rotated_polygon <- spatialEco::rotate.polygon(polygon_sf, angle = angle, anchor = "center")
#
# # Extract the geometry from the rotated sf object
# return(st_geometry(rotated_polygon)[[1]])
# })
#
# # Combine the results into a single sf object
# buffered_polygons_sf <- sf::st_sfc(buffered_polygons_list, crs = sf::st_crs(seeded_points)) |> st_as_sf()
#
# # Intersect the plots with the reef polygons, calculate overlap
# intersections <- st_intersection(buffered_polygons_sf, st_union(geomorphic_seeding)) %>%
# mutate(intersection_area = st_area(.)) %>%
# mutate(percentage_overlap = (as.numeric(intersection_area) / 10000) * 100)
#
# return(intersections)
# })
#
# # Combine all results into one sf object
# final_results <- do.call(rbind, results_list)
#
# # Reset the plan to sequential processing after completion (optional)
# plan(sequential)
#
# toc()
#
# # Final results
# final_results
library(future.apply)
library(sf)
library(spatialEco)
tic()
# Set up parallel processing plan
plan(multisession, workers = availableCores())
# Use future_lapply to create the buffered polygons with random rotations
buffered_polygons_list_rotated <- future_lapply(seq_len(nrow(seeded_points)), function(i) {
# Extract the coordinates of the current point
x <- sf::st_coordinates(seeded_points)[i, 1]
y <- sf::st_coordinates(seeded_points)[i, 2]
# Set parameters for the rectangle around the point
x_min <- x - (width / 2)
x_max <- x + (width / 2)
y_min <- y - (length / 2)
y_max <- y + (length / 2)
# Create the rectangular polygon
polygon <- sf::st_polygon(list(rbind(
c(x_min, y_min),
c(x_min, y_max),
c(x_max, y_max),
c(x_max, y_min),
c(x_min, y_min)
)))
# Convert to sf object for rotation
polygon_sf <- st_sf(geometry = st_sfc(polygon, crs = st_crs(seeded_points)))
# Select a random rotation angle from the sequence
angle <- sample(seq(0, 355, 5), 1)
# Rotate the polygon by the random angle using spatialEco::rotate.polygon
rotated_polygon <- spatialEco::rotate.polygon(polygon_sf, angle = angle, anchor = "center") |> mutate(id=i)
# Extract the geometry from the rotated sf object
return(rotated_polygon)
})
buffered_polygons_rotated <- do.call(rbind,buffered_polygons_list_rotated) |> st_set_crs(20353)
buffered_polygons_rotated_intersect <- st_intersection(buffered_polygons_rotated, st_union(geomorphic_seeding)) %>%
mutate(intersection_area = st_area(.)) |> as.data.frame() |> select(-geometry)
intersections_rotated <- buffered_polygons_rotated |>
left_join(buffered_polygons_rotated_intersect, by="id") |>
mutate(percentage_overlap = (as.numeric(intersection_area) / 10000) * 100)
# Reset the plan to sequential processing after completion (optional)
plan(sequential)
toc()
## 45.9 sec elapsed
tmap_mode("plot")
tm_basemap("Esri.WorldImagery") +
tm_shape(intersections_rotated |> filter(!percentage_overlap>98),
bbox=inset) +
tm_polygons(fill=NA,
lwd=0.5,
col_alpha=0.5,
col="darkred") +
tm_shape(intersections_rotated |> filter(percentage_overlap>98),
bbox=inset) +
tm_polygons(fill=NA,
lwd=1,
col="white")